home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Fritz: All Fritz
/
All Fritz.zip
/
All Fritz
/
FILES
/
PROGMISC
/
PCSSP.LZH
/
PC-SSP.ZIP
/
MATINSL.ZIP
/
FACTR.FOR
< prev
next >
Wrap
Text File
|
1985-11-29
|
3KB
|
128 lines
C
C ..................................................................
C
C SUBROUTINE FACTR
C
C PURPOSE
C FACTORIZATION OF THE MATRIX A INTO A PRODUCT OF A LOWER
C TRIANGULAR MATRIX L AND AN UPPER TRIANGULAR MATRIX U. L HAS
C UNIT DIAGONAL WHICH IS NOT STORED.
C
C USAGE
C CALL FACTR(A,PER,N,IA,IER)
C
C DESCRIPTION OF PARAMETERS
C A MATRIX A
C PER ONE DIMENSIONAL ARRAY WHERE PERMUTATIONS OF ROWS OF
C THE MATRIX ARE STORED
C DIMENSION OF PER MUST BE GREATER THAN OR EQUAL TO N
C N ORDER OF THE MATRIX A
C IA SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY A
C IN THE CALLING PROGRAM WHEN THE MATRIX IS IN DOUBLE
C SUBSCRIPTED DATA STORAGE MODE. IA=N WHEN THE MATRIX
C IS IN SSP VECTOR STORAGE MODE.
C IER ERROR INDICATOR WHICH IS ZERO IF THERE IS NO ERROR,
C AND IS THREE IF THE PROCEDURE FAILS.
C
C REMARKS
C THE ORIGINAL MATRIX, A,IS REPLACED BY THE TRIANGULAR FACTORS
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C SUCCESSIVE COMPUTATION OF THE COLUMNS OF L AND THE
C CORRESPONDING ROWS OF U.
C
C REFERENCES
C J. H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM -
C CLARENDON PRESS, OXFORD, 1965. H. J. BOWDLER, R. S. MARTIN,
C G. PETERS, AND J. H. WILKINSON - 'SOLUTION OF REAL AND
C COMPLEX SYSTEMS OF LINEAR EQUATIONS', NUMERISCHE MATHEMATIK,
C VOL. 8, NO. 3, 1966, P. 217-234.
C
C ..................................................................
C
SUBROUTINE FACTR(A,PER,N,IA,IER)
DIMENSION A(1),PER(1)
DOUBLE PRECISION DP
C
C COMPUTATION OF WEIGHTS FOR EQUILIBRATION
C
DO 20 I=1,N
X=0.
IJ=I
DO 10 J=1,N
IF (ABS(A(IJ))-X)10,10,5
5 X=ABS(A(IJ))
10 IJ=IJ+IA
IF (X) 110,110,20
20 PER(I)=1./X
I0=0
DO 100 I=1,N
IM1=I-1
IP1=I+1
IPIVOT=I
X=0.
C
C COMPUTATION OF THE ITH COLUMN OF L
C
DO 50 K=I,N
KI=I0+K
DP=A(KI)
IF (I-1) 110,40,25
25 KJ=K
DO 30 J=1,IM1
IJ=I0+J
DP=DP-1.D0*A(KJ)*A(IJ)
30 KJ=KJ+IA
A(KI)=DP
C
C SEARCH FOR EQUILIBRATED PIVOT
C
40 IF (X-DABS(DP)*PER(K))45,50,50
45 IPIVOT=K
X=DABS(DP)*PER(K)
50 CONTINUE
IF (X)110,110,55
C
C PERMUTATION OF ROWS IF REQUIRED
C
55 IF (IPIVOT-I) 110,70,57
57 KI=IPIVOT
IJ=I
DO 60 J=1,N
X=A(IJ)
A(IJ)=A(KI)
A(KI)=X
KI=KI+IA
60 IJ=IJ+IA
PER(IPIVOT)=PER(I)
70 PER(I)=IPIVOT
IF (I-N) 72,100,100
72 IJ=I0+I
X=A(IJ)
C
C COMPUTATION OF THE ITH ROW OF U
C
K0=I0+IA
DO 90 K=IP1,N
KI=I0+K
A(KI)=A(KI)/X
IF (I-1)110,90,75
75 IJ=I
KI=K0+I
DP=A(KI)
DO 80 J=1,IM1
KJ=K0+J
DP=DP-1.D0*A(IJ)*A(KJ)
80 IJ=IJ+IA
A(KI)=DP
90 K0=K0+IA
100 I0=I0+IA
IER=0
RETURN
110 IER=3
RETURN
END